Energetics of hydrogen/lithium complexes in silicon analyzed using the Maxwell construction 
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We have studied hydrogen/lithium complexes in crystalline silicon using density-functional-theory methods 
and the ab initio random structure searching (AIRSS) method for predicting structures. A method based on 
the Maxwell construction and convex hull diagrams is introduced which gives a graphical representation of 
the relative stabilities of point defects in a crystal and enables visualization of the changes in stability when 
the chemical potentials are altered. We have used this approach to study lithium and hydrogen impurities in 
silicon, which models aspects of the anode material in the recently-suggested lithium-ion batteries. We show 
that hydrogen may play a role in these anodes, finding that hydrogen atoms bind to three-atom lithium clusters 
in silicon, forming stable {H,3Li} and {2H,3Li} complexes, while the {H,2Li} complex is almost stable. 
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I. INTRODUCTION 

The properties of a material may be substantially affected 
by the presence of impurity atoms. 1 One of the fundamental 
properties is the relative energies of different configurations 
of the impurities within the material. Several different species 
of atom may be present, and the impurity atoms may have 
strong interactions with each other and with the host mate- 
rial. In this paper we study the interactions between lithium 
and hydrogen impurities in crystalline silicon, which is a sys- 
tem of interest for lithium-ion battery technologies. Lithium- 
ion batteries are widely used in portable electronics and are 
now being employed in the next generation of hybrid and all- 
electric vehicles. 2 Recently there has been interest in using 
silicon anodes, which have a very high volumetric and gravi- 
metric capacity. 3 " 5 Crystalline silicon is generally used in the 
battery and the first charge-discharge cycle involves the con- 
version of crystalline silicon to an amorphous phase. 

Point defects may form a variety of complexes containing 
several atoms of more than one atomic species. An interesting 
feature of our approach is the use of a graphical representation 
of the relative stabilities of the defects based on the Maxwell 
construction 6 - 7 and the corresponding convex hull. This is a 
standard procedure used in materials science (and elsewhere) 
for studying phase separation, which we have applied to de- 
fects in crystals. The term "convex hull" refers to the fact that 
the second law of thermodynamics requires that the energy per 
particle (or free energy as appropriate) must be a convex func- 
tion of the relative concentrations of the particles. Appropriate 
chemical potentials for the different atomic species can nor- 
mally be estimated for a particular set of external conditions. 
We are also interested in the changes in the relative stabilities 
of defects when the chemical potentials of the atomic species 
are altered. Charged defects are also of interest in general, and 
these can be included by introducing a chemical potential for 



electrons, although we will not show results for this case here. 

This article is divided into the following sections, in Sect. [II] 
the modified Maxwell construction is introduced to visualize 
the energetics of the phase separation of defect complexes, 
and the example of N/O complexes in silicon is used to vali- 
date the method. In Sect.|III]the method is applied to the prob- 
lem of H/Li defect complexes in silicon. Finally Sect. lIVI con- 
cludes the work. 



II. CONVEX HULL FOR DEFECTS 

In this article we denote a defect complex by listing its con- 
stituent atoms between braces, for example, the {2H,3Li} de- 
fect complex contains two hydrogen and three lithium atoms. 
For simplicity we begin by considering a single impurity 
species a which forms a single impurity complex in a large 
sample of perfect bulk crystal. The formation energy of a 
complex {n a a} containing n a impurity atoms is given by 



If = E D - n a H a - E h 



(1) 



where E^, is the energy of the system including the defect 
complex, fi a is the chemical potential of the impurity species, 
and En is the energy of the host crystal without defects. The 
energy of the host crystal can be written as 



Eh = Y< n pVp> 



(2) 



where there are n$ atoms of the host crystal atoms of each 
species /3 which have chemical potentials The complex 
with the lowest value of Ef/n a is the most stable and ones 
with higher energies are metastable. 

Eq. (Q~|i can readily be extended to a complex containing a 
larger number of impurity species, and the formation energy 
per impurity atom is then 
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The convex hull diagram is constructed by plotting the for- 
mation energy per impurity atom of each complex against the 
fractional concentration C, of impurity species i, where 



If the complexes contain only two atomic species we can sim- 
ply plot £p a (C,) on a graph as shown, for example, in Fig. 
Q] The solid red and blue lines in Fig. Q]show convex hulls 
which consist of straight-line segments known as "tie lines". 
A higher-dimensional hull is required if there are more than 
two impurity species. The case of three atomic species can be 
represented on a "triangular plot" as shown, for example, in 
Ref.[H Alternatively one can consider two-dimensional slices 
of the multi-dimensional hull. Electrons can also be included 
as a species, so that defect complexes in different charge states 
can be described. 

When more than one impurity species is present it is possi- 
ble for "phase separation" or "disproportionation" of defect 
complexes to occur. For example, the lowest energy state 
when equal numbers of impurity species A and B are present 
might consist of equal numbers of the defect complexes 
{A,2B} and {A}, rather than the defect complex {2A,2B}. 
We will assume that the defect complexes do not interact with 
one another (except by a "chemical reaction", forming new 
complexes which can be included within the theory) and they 
are expected to be distributed randomly over the host crystal 
rather than forming spatially distinct phases as in the standard 
application to "phase separation". We may include within E va 
the contributions to the free energy from the vibrational en- 
tropy and the configurational entropy due to the different de- 
generate orientations of a complex at a lattice site. However, 
it is not possible to include the configurational entropy due 
to arranging defects at different lattice sites within £p a , since 
it depends on the concentrations of the impurity complexes. 
Hence the Maxwell construction is exact only at zero temper- 
ature, and a full statistical mechanical description is required 
at finite temperatures. The standard plot of the defect forma- 
tion energy against the chemical potential found in the defect 
literature (see, for example, Fig. 2 of Ref. 9), suffers from the 
same limitation. 



A. Application to N/O complexes in silicon 

We use N/O complexes in silicon to illustrate our approach 
and outline the advantages of the method. Fig.[T]shows a con- 
vex hull containing N/O complexes from a previous study. 10 
The hull shown in the upper panel of Fig. Q] is constructed as 
follows. We choose the chemical potentials for the N and O 
atoms to be the same as in our previous study. 10 The chemical 
potential for N is obtained from the energy of the {2N} defect 
as 

Al N ({2N}) = ±(£ D ({2N})-£ H ), (5) 
while for O it is obtained from the energy of the {O} defect as 
Ho({0})=E»({0})-E H . (6) 



These energies are plotted on the diagram at the points 
(£pa=0,CN=0) and (0,1). The solid-red convex hull is then 
drawn as a series of straight-line segments between (0,0) and 
(0, 1) passing through the squares representing the energies per 
impurity atom £p a in such a way that none of the defects lie 
below the hull. Only the complexes which lie on the hull are 
stable. The {O}, {N,30}, {N,20}, {2N,20} and {2N,Oj 
defects are all stable since they lie on the tie line of the con- 
vex hull, but the {N,20} complex is not stable and the energy 
would be lowered if these defects disproportionated such that 



4{N,20} ^2{N,30} + {2N,20}. (7) 

The dashed blue line in the upper panel of Fig.[T]indicates 
how the situation would be modified if the chemical poten- 
tial for O is changed to the energy of the much more stable 
{40} complex. The modified complex hull is redrawn as the 
solid-blue line in the lower panel of Fig. Q] It is clear from 
either the upper or lower panels that the {N,30} defect is no 
longer stable under these circumstances and the energy would 
be lowered if they disproportionated such that 

2{N,30}^{40} + {2N,20}. (8) 

These results are entirely consistent with the analysis reported 
in Ref. 10. The graphical representation is useful because it 
shows which defects are stable for a particular choice of chem- 
ical potentials, and it also allows an easy visualization of the 
effects of altering the chemical potentials. The convex hull 
also shows the amount by which the energy is lowered by the 
disproportionation of defect complexes. 



B. Native defects 

Our approach is easily extended to account for a crys- 
tal containing n v native vacancies V which form com- 
plexes {n a (X,npl3,n v V} or «/ interstitials / forming complexes 
{n a (X,npl5,niI}. The native defects may be included within 
this framework, not as new "impurities species", but as a 
change in the definition of the host material. The chemical 
potential of an atom of the host material added at a surface of 
the crystal is equal to the cohesive energy of an atom in the 
perfect bulk crystal since and, as long as the surface is large, 
the surface energy is unchanged by adding an atom at a "kink 
site" on a surface terrace. Hence in the case where the vacancy 
is absent prior to the addition of impurities, fip is the chemical 
potential of each host species /3 in the perfect crystal. If the 
vacancy is present prior to the addition of impurities, jip is the 
chemical potential of host species j3 in a crystal containing the 
vacancy. Note that it is not possible to describe more than one 
type of vacancy in this approach unless each type is treated as 
a separate species. The same applies for interstitials. 
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FIG. 1: (Color online) Maxwell construction for N/O complexes 
in silicon. The upper panel shows the formation energy per impu- 
rity atom (£p a ) for various complexes using the {2N} and {0} de- 
fects for the chemical potentials. The stable complexes, {N,30}, 
{2N,20} and {2N,0} are found on the convex hull (solid red). The 
blue (dashed) line is a tie-line between the values of £ pa for the {40} 
and {2N} complexes. If the O chemical potential is changed to that 
of the more stable {40} complex, the stable complexes become those 
which both lie on the solid-red convex hull and also lie below the blue 
dashed line. This construction shows the stable complexes when the 
O chemical potential is changed to the energy of {40}, since only 
the complexes with formation energies below the tie line ({2N,20} 
and {2N,0}) are stable. This is confirmed by replotting the energies 
(blue diamonds) in the lower panel using {40} as the chemical po- 
tential. Only two structures lie on the new convex hull (solid blue 
line). 



III. STUDY OF H/LI COMPLEXES IN SILICON 

A. H and Li in silicon 

Hydrogen is a common impurity in silicon which binds to a 
variety of defect complexes, including those containing oxy- 
gen and nitrogen^ It is difficult to determine the concen- 
tration of hydrogen in silicon samples although it is likely 
to be present at virtually every stage of manufacture and it 
may be incorporated at levels of up to about 1 20 cm" 1 J^- The 
role of hydrogen impurities in semiconductors has been re- 
viewed by Estreicher. 12 In this article we present results for 
H/Li impurities in bulk silicon. Both amorphou a 13 ' 14 and 
crystalline 15 phases of the Li-Si system have been studied us- 
ing density-functional-theory (DFT) methods, giving voltages 
which agree with experiment to within 0. 1 V. Wan et alJ& 
studied lithium defects in silicon using DFT and found the 



most stable position for a single lithium impurity to be at the 
tetrahedral (T c i) site, which is preferred to substitution on a sil- 
icon site. Kim et ah xl performed a detailed study of the elec- 
tronic structure of a lithium atom in bulk silicon, concluding 
that the presence of lithium weakens the nearby Si-Si bonds. 

B. DFT calculations 

We have predicted structures of lithium impurities in sil- 
icon using DFT and the ab initio random structure search- 
ing method (AIRSS). 18 ' 19 In this method randomly gener- 
ated structures are relaxed to a minimum in the energy. This 
approach has been successful in predicting the structures of 
point defects in semiconductors 10,20 and ceramics, 2 ' as well 
as high-pressure phases of materials such as solid silane, 18 
hydrogen, 22 and lithium. 23 Each AIRSS run was performed 
at a fixed stoichiometry. Searches were performed at several 
stoichiometries and the results were analyzed using the con- 
vex hull construction. 

We used simulation cells containing 1, 2, 3 or 4 lithium 
atoms, 1, 2, 3 or 4 lithium atoms plus 1 hydrogen, atom and 
1, 2 or 3 lithium plus 2 hydrogen atoms. The initial structures 
were chosen in a similar fashion to those in our earlier study 
of H/N/O defects. 10 A small hole was made in the host crystal 
by removing a silicon atom from the unit cell and placing the 
appropriate number of lithium atoms and a silicon atom ran- 
domly within the hole. The configurations were then relaxed 
using the DFT forces. Some of the starting configurations re- 
laxed to structure in which two or more separate defects were 
present, and these structures were discarded. 

The plane-wave basis-set DFT code CASTEp2i was used 
to calculate ground-state total energies and optimize the ge- 
ometries of the candidate structures. We used "on-the-fly" 
ultrasoft pseudopotentials 25 and the PBE 26 parameterization 
of the generalized gradient approximation to the exchange- 
correlation functional. The initial searches were performed 
with simulation cells containing 32 silicon atoms plus im- 
purity atoms, sampling the Brillouin zone using the 2 x 
2x2 multi-fc-point generalization 20 ' 27 ' 28 of the Baldereschi 
scheme. 29 The basis set contained all plane-waves with en- 
ergies up to 300 eV. 

The lowest-energy structures from each search were then 
further relaxed at a higher level of accuracy. We used simula- 
tion cells containing 256 silicon atoms and a basis set with 
plane-waves up to 500 eV, a harder Li pseudopotential and 
a standard 2x2x2 Monkhorst-Pack grid of ^-points. The 
formation energies were then calculated using these struc- 
tures but with a larger 4x4x4 Monkhorst-Pack grid of k- 
points. We also performed some calculations in which we 
allowed spin polarization to occur, but no significant effects 
were found. 



C. Results for H/Li in silicon 

We model the initial stages of the incorporation of Li atoms 
into bulk silicon as the formation of defect complexes with 
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Defect 


Ef (eV) 


{3Li} 


0.422 


{H,3Li}* 


0.103 


{H,3Li} 


-0.139 


{H,2Li}** 


0.321 


{H,2Li}* 


0.301 


{H,2Li} 


0.101 


{2H,3Li}** 


0.114 


{2H,3Li}* 


-0.010 


{2H,3Li} 


-0.239 


{H,Li}** 


0.474 


{H,Ll| 




{H,Li} 


0.401 


{2H,2Li} 


-0.028 


{2H,Li}* 


0.413 


{2H,Li} 


-0.076 



TABLE I: Formation energies, Ef, of various H/Li complexes in sil- 
icon. The chemical potential for H is obtained from the energy of 
an H2 molecule at the site of bulk silicon and the chemical po- 
tential for Li is obtained from the energy of a Li atom at the Tj site. 
Single asterisks (*) and double asterisks (**) denote the first and sec- 
ond metastable defects, respectively, for a given defect stoichiometry. 
The same information is represented graphically in Fig.|2] 




FIG. 3: (Color online) Two views of the {H,3Li} defect in bulk sili- 
con. Silicon atoms are shown in yellow (light grey), lithium in purple 
(dark grey) and hydrogen in white. This defect is predicted to be sta- 
ble for < Ch < 0.4. The lithium atoms form a three-membered 
ring in a void created by breaking a Si-Si bond and a hydrogen atom 
bonds to the displaced silicon atom. 
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FIG. 2: (Color online) Convex hull for Li/H complexes in silicon. 
Formation energy per impurity atom plotted against the fractional 
concentration Ch of hydrogen. The plus signs (+) represent defect 
complexes. The chemical potentials for Li ({Li}) and H ({H2}) are 
plotted at (0,0) and (0,1), respectively. Metastable complexes of a 
given stoichiometry are denoted by asterisks (*). The formation en- 
ergies of defect complexes are indicated by blue (dark grey) crosses 
and the predicted stable defects by red (mid-grey) crosses. The red 
(mid-grey) lines show the formation energy per impurity atom of the 
(mixtures of) stable defects present as a function of Ch- 



FIG. 4: (Color online) The {2H,3Li} complex in bulk silicon. Sili- 
con atoms are shown in yellow (light grey), lithium in purple (dark 
grey) and hydrogen in white. This defect is predicted to be stable 
for 0.25 < Ch < L0. It is similar to the {H,3Li} defect shown in 
Fig.[U except that the additional hydrogen atom terminates the dan- 
gling bond on the displaced silicon atoms. The lithium atoms form 
a three-membered ring in the void created by breaking a Si-Si bond. 
A hydrogen atom sits between these two silicon atoms. 



varying concentrations of H and Li atoms. We assume that 
the Li atoms are supplied from bulk body-centered-cubic Li 
metal. The {Li} defect in silicon has a formation energy of 
0.21 eV relative to bulk Li in accordance with the observed 
overpotential during Li insertion reactions in bulk silicon. 32 
We assume that the Li atoms initially form {Li} defects in 
the bulk silicon and we investigate whether chemical reactions 
between Li and H atoms in bulk silicon are energetically fa- 
vorable. We choose the {Li} defect to define our chemical 
potential because it is the lowest energy defect that we have 
found in silicon and because a lithium defect of the same {T c f) 
symmetry has been observed in electron paramagnetic reso- 
nance experiments in bulk silicon^ A metastable {3Li} com- 
plex with Ci v symmetry was also found with formation energy 
relative to {Li} of 0. 14 eV per atom, but this is too high in en- 
ergy to appear on Fig. [2] The hydrogen chemical potential 
was obtained as in previous studies 10 - 20 from a H2 molecule 
at the T c i site in silicon. A more complete model would require 
the inclusion of other defects, such as vacancies and vacancy 
complexes. 

Single-atom Li defects are favored when only Li impurities 
are present, but the convex hull plotted in Fig.|2]shows that the 
addition of hydrogen causes the lithium atoms to cluster. For 
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< Ch < 0.25 the {H,3Li} defect forms and any additional 
lithium atoms occupy T c / sites. The structure of the {H,3Li} 
defect is pictured in Fig.[3] Both the {H,3Li} and {2H,3Li} 
complexes are present in the range 0.25 < Ch < 0.40. The 
{2H,3Li} defect is present when Ch > 0.40, and any excess 
hydrogen forms H2 molecules at the T c / site. The {H,3Li}, 
{2H,3Li} and {2H,Li} defects resemble a Chadi defect 31 
with the lithium atoms close by. The metastable com- 
plex found by Chadi 3 1 has a formation energy relative to a H2 
molecule at the Tj site of 0. 1 eV per H atom. 20 The {2H,Li} 
defect could be present when thermal effects are taken into 
account, due to its close proximity to the tie line at Ch = 0.66. 

IV. CONCLUSIONS 

We have used convex hull diagrams based on the Maxwell 
construction to visualize the energetics of defect complexes. 
The graphical representation shows the stable defects for a 
particular choice of chemical potentials and the relative sta- 
bilities of the defect complexes. It allows visualization of the 
changes in the stability of the defects when the chemical po- 
tentials are changed. The approach was illustrated using data 
for N/O complexes in silicon from a previous study.— The 
method can easily be extended to crystals containing charged 
species, native vacancies or interstitials, which is an area we 
are currently investigating. 



We also studied H/Li complexes in silicon, which may be 
relevant to understanding how the first few Li atoms enter 
crystalline silicon. This extended the previous analysis of 
lithium in silicon anodes to include its interaction with the hy- 
drogen that is introduced into the silicon during manufacture. 
When only lithium impurity atoms are present the most stable 
defect is {Li}, but the presence of hydrogen atoms causes the 
lithium atoms to form H/Li clusters containing three lithium 
atoms. We found that the {H,3Li} and {2H,3Li} complexes 
are stable and that the {2H,Li} complex is nearly stable. Hy- 
drogen present in the bulk favors bonding covalently to self- 
interstitial defects, 20 but once the silicon defects are saturated, 
hydrogen will bind to lithium complexes, further stabilizing 
them. The combination of AIRSS and the Maxwell construc- 
tion and convex hull diagrams is a powerful tool for analyzing 
the energetics of point defects in materials. 
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